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Abstract. In this paper, we provide a simple framework to derive and analyse several classes of effective one- 
step methods. The framework consists in the discretization of a local Fourier expansion of the continuous problem. 
Different choices of the basis lead to different classes of methods, even though we shall here consider only the case 
of an orthonormal polynomial basis, from which a large subclass of Runge-Kutta methods can be derived. The 
obtained results are then applied to prove, in a simplified way, the order and stability properties of Hamiltonian 
BVMs (HBVMs), a recently introduced class of energy preserving methods for canonical Hamiltonian systems (see 
[l] and references therein). A few numerical tests with such methods are also included, in order to confirm the 
effectiveness of the methods. 
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1. Introduction. 

Though I have not always been able to make simple 
a difficult thing, I never made difficult a simple one. 
F. G. Tricomi 

One-step methods are widely used in the numerical solution of initial value problems for ordinary 
differential equations which, without loss of generality, we shall assume to be in the form: 

y'(t) = f(y(t)), t€[t ,t +T], y(t Q ) = y Q eR m . (1.1) 

In particular, we consider a very general class of effective one-step methods that can be led back 
to a local Fourier expansion of the continuous problem over the interval [to, to + h], where h is the 
considered stepsize. In general, different choices of the basis result in different classes of methods, 
for which, however, the analysis turns out to be remarkably simple. Though the arguments can be 
extended to a general choice of the basis, wc consider here only the case of a polynomial basis, from 
which one obtains a large subclass of Runge-Kutta methods. Usually, the order properties of such 
methods are studied through the classical theory of Butcher on rooted trees (see, e.g., [H Th. 2.13 
on p. 153]), almost always resorting to the so called simplifying assumptions (see, e.g., Th. 7.4 on 
p. 208]). Nonetheless, such analysis turns out to be greatly simplified for the methods derived in the 
new framework, which is introduced in Section [2] Similar arguments apply to the linear stability 
analysis of the methods, here easily discussed through the Lyapunov method. Then, we apply the 
same procedure to the case where (|1.1[) is a canonical Hamiltonian problem, i.e., a problem in the 
form 

^ = JVH(y), J = f J 7 ™ ) , y(t ) = yo & R 2 "\ (1.2) 
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where H(y) is a smooth scalar function, thus obtaining, in Section [3J an alternative derivation of 
the recently introduced class of energy preserving methods called Hamiltonian BVMs (HBVMs, 
see [U [21 [3] and references therein). A few numerical examples concerning such methods are then 
provided in Section [4l in order to make evident their potentialities. Some concluding remarks are 
then given in Section [5] 

2. Local Fourier expansion of ODEs. Let us consider problem (jl.ll) restricted to the interval 
[t ,t + h]: 

y' = f(y), te[to,to + h], y(t ) = yo. (2.1) 

In order to make the arguments as simple as possible, we shall hereafter assume / to be analytical. 
Then, let us fix an orthonormal basis {Pj}°^ Q over the interval [0,1], even though different bases 
and/or reference intervals could be in principle considered. In particular, hereafter we shall consider 
a polynomial basis: i.e., the shifted Legendre polynomials over the interval [0, 1], scaled in order to 
be orthonormal. Consequently, 

/ P l (x)Pj (x) dx = 8 lj , deg Pj = j, Vi, j > 0, 
Jo 

where 5ij is the Kroneckcr symbol. We can then rewrite (|2.1I) by expanding the right-hand side: 
oo „i 
y'(t + ch) = Y^Pj(chj(y), ce[0,l]; 7j (y) = / P,(r)/(j/(t + rh)) dr. (2.2) 

The basic idea (first sketched in [5]) is now that of truncating the series iafter r terms, which turns 
([2~2]) into 

r-1 „1 

w% + ch)=J2P j (c)>y j (w), c€[0,l]; lM = Pj(r)f(w(t +rh))dr. (2.3) 

By imposing the initial condition, one then obtains 

uj(to + ch)=yo + hy2lj(w) P j (x)dx, ce[0,l]. (2.4) 

Obviously, us is a polynomial of degree at most r. The following question then naturally arises: 
"how close are y(to + h) and u)(to + h) ?" The answer is readily obtained, by using the following 
preliminary result. 

Lemma 2.1. Let g : [0, h] — > W n be a suitably regular function. Then J" 1 Pj(T)g(Th) dr = 0(h j ). 
Proof. Assume, for sake of simplicity, 

n=0 

to be the Taylor expansion of g. Then, for all j > 0, 

rl 00 n (n)(r\\ rl 

/ Pj ( T )g( Th )dT = J2 L -± 1 h n / P j (T)T n dT = 0(V), 
Jo n=0 nl Jo 



since Pj is orthogonal to polynomials of degree n < j. □ 

As a consequence, one has that (see (|2.3p ) = 0(/i J ). Moreover, for any given t £ [t ,t + h], 

we denote by y(s,i,y) the solution of (|2.1[) - (|2.2p at time s and with initial condition y(t) = y. 
Similarly, we denote by 

$(s,t,y) = ^y(s,t,y), (2.5) 

also recalling the following standard result from the theory of ODEs: 

JLy(«, t, y) = -*(«, t, (2.6) 
at 

We can now state the following result, for which we provide a more direct proof, with respect 
to that given in [5). Such proof is essentially based on that of [HI Theorem6.5.1 on pp. 165-166]. 

Theorem 2.2. Let y(t Q + ch) and uj(t + ch), c £ [0,1], be the solutions of h2.2}) and h2.S\) , 
respectively. Then, y(tg + h) — ui(to + h) = 0{h 2r+1 ). 

Proof. By virtue of Lemma 12.11 and (|2.5[) - (|2.6p . one has: 

y(t + h) -u(t + h) = y(t + h, t , y ) - y(to + h,t + h, u(t + h)) 

fto+h d Pt + h / Q Q X 

— y(t + h,T,u(r))dT = J \^—y(t + h,T,u(T)) + —y(t Q + h,T 1 uj(T))uj'(T)JdT 
h [ $(t + h,t + ch,uj(t + ch)) {-f(u(t + ch)) + oj'{t + ch)) dc 



-h [ <P(t + h,t + ch,uj{t + chj) ( V 7j (o;)P J (c) ] dc 
J ° V=r ^ J 

OO / n\ \ OO 

' h ^\J Pj^Mto + ^tQ + ch^ito + ch^dcJj^Lj) = h^O{W)0{h J ) = 0(h 2r+1 ). 



J=r " j=r 



The previous result reveals the extent to which the polynomial ui(t), solution of p.3p . approxi- 
mates the solution y(t) of the original problem (|2.1j) on the time interval [to, to + h]. Obviously, the 
value w(to +h) may serve as the initial condition for a new IVP in the form (|2.3j) approximating y(t) 
on the time interval [to + h, to + 2h]. In general, setting U = to + ih, i = 0, 1, . . . , and assuming that 
an approximation uj(t) is available on the interval [ti-2,U-i], one can extend the approximation to 
the interval [tj_i,tj] by solving the IVP 



r-l „i 

'(U- 1 +ch)=J2P J (c) P,(r)/(o;(t i _ 1 +c/i))dr, ce[0,l], 



(2.7) 



the initial value cj(£j_i) having been computed at the preceding step. The approximation to y(t) is 
thus extended on an arbitrary interval [to, to + Nh), and the function ui(t) is a continuous piecewise 
polynomial. As a direct consequence of Theorem 12. 2[ we obtain the following result. 



Corollary 1. Let T = Nh, where h > and N is an integer. Then, the approximation to 
the solution of problem by means of (|2.7p at the grid-points ti = £^_i + h, % = 1, . . . ,N, with 

w (^o) = Ho, is 0{h 2r ) accurate. 
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We now want to compare the asymptotic behavior of u(t) and y(t) on the infinite length interval 
[t , +00) in the case where / is linear or defines a canonical Hamiltonian problem. To this end we 
introduce the the infinite sequence {oJi\ = {w(ti)}. 

Remark 1. Though in general, the sequence {w^} cannot be formally regarded as the outcome 
of a numerical method, under special situations, this can be the case. For example, when f is a 
polynomial, the integrals in (|2.3[) may be explicitly determined and the IVP in ()2.3|) is evidently 
equivalent to a nonlinear system having as unknowns the coefficients of the polynomial u expanded 
along a given basis (for example, the polynomyal uj may be computed by means of the method of 
undetermined coefficients). This issue, as well as details about how to manage the integrals in 
the event that the integrands do not admit an analytical primitive function in closed form, will be 
thoroughly faced in Section \3\ 

2.1. Linear stability analysis. For the linear stability analysis, problem (|2.1[) becomes the 
celebrated test equation 



y' = Ay, R(A) < 0. 



(2.8) 



By setting 



X = a + i(3, y = xi+ix 2 , 

with i the imaginary unit, problem 

-' - Ax 



(xi,x 2 ) 1 



A 



a ~P 
P a 



can be rewritten as 
x' = Ax, t £ [to, t + h], x(to) given. 
Consequently, the corresponding truncated problem (|2.3[) becomes 

r-1 .1 

u'(to + ch) = ^VPj(c) / Pj(T)VV(oj(to + Th))dT, ce[0,l], 



3=0 



where 



is a Lyapunov function for 



V(x) 



- T 
-X X 



(2.9) 



(2.10) 



(2.11) 



From (|2.10|) - (|2.1ip one readily obtains 



AV(w(t )) = V(u(t + h)) -V(w(to)) = hi VV(u{t + rh)) 1 w'(t + rh) dr 

r-1 r /.l 1 J T 

P 3 (T)\7V{0j(t +Th))dT A I 



r-1 r r l 

,._n Uo 



3=0 
r-1 



Pj(T)VV(uj(to+Th))dT 



r-1 „i 

3=0 J ° 



Pj(T)w(t +Th) dr 



Last equality follows by taking the symmetric part of A. We observe that 



U3 ^ 



E 

3=0 



Pj(T)uj(t + Th) dr 
4 



> 0, 



since, conversely, this would imply co(to + ch) = p ■ P r (c) for a suitable p ^ and, therefore (from 
(I2.10[) ), P' r = which is clearly not true. Thus for a generic y ^ 0, 

Ay(w(t )) < ^ »(A) < and AV{u(t )) = K(A) = 0. 

Again, the above computation can be extended to any interval [ti—i, U] and, from the discrete version 
of the Lyapunov theorem (see, e.g., [HI Th. 4.8.3 on p. 108]) , we have that the sequence Wj tends to 
zero if and only if 5R(A) < 0, while it remains bounded whenever Ji(A) = 0, whatever is the stepsizc 
h > used. The following result is thus proved. 

Theorem 2.3. The continuous solution y(t) of (|2.8|) and its discrete approximation Ui have 
the same stability properies, for any choice of the stepsize h > 0. 

2.2. The Hamiltonian case. In the case where problem ()1.1|) is Hamiltonian, i.e., (|1.2j) . the 

approximation provided by the polynomial uj in (|2.3p - (|2.4[) inherits a very important property of the 
continuous problem, i.e., energy conservation. Indeed, it is very well known that for the continuous 
solution one has, by virtue of (|1.2[) , 

^H{y{t)) = VH(y(t)) T y'(t) = VH(y(t)) T JVH(y(t)) = 0, 

due to the fact that matrix J is skew-symmetric. Consequently, H(y(t)) = H(y ) for all t. For the 
truncated Fourier problem, the following result holds true. 

Theorem 2.4. H(ui(t + h)) = H(u(t )) = H(y ). 

Proof. From (|2.3[) . considering that f(uj) = J\7H(uj) and J T J = I, one obtains: 
H(u(to + h))-H(yo) = 

pi pi r-1 

= k VH(uj(t a +Th)) T Uj'(t +Th)dT = hi VH(uj(t +Th)) T ^P ] {T) lj {uj)dT 

Jo Jo j=0 

= h ll{l VHMto + TZO^todr) 7» = hJ^^iufJ^iu) = 0, 

since J is skew-symmetric. □ 

3. Discretization. Clearly, the integrals in ()2.3|) . if not directly computable, need to be nu- 
merically approximated. This can be done by introducing a quadrature formula based at k > r 
abscissae {c^}, thus obtaining an approximation to (|2.3[) : 

r-1 fe 

u\t + ch)=^P j (c)^btP j (c l )f(u(to + c e h)), ce[0,l]. (3.1) 

3=0 e=i 

where the {bt} are the quadrature weights, and u is the resulting polynomial, of degree at most r, 
approximating ui. It can be obtained by solving a discrete problem in the form: 

r-1 k 

u'(t + c i h) = J2P j (c i )Y, b ePj(ce)f(u(to + c £ h)), i = l,...,k. (3.2) 

3=0 1=1 
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Let q be the order of the formula, i.e., let it be exact for polynomials of degree less than q (we 
observe that q > k > r). Clearly, since we assume / to be analytical, choosing k large enough, along 
with a suitable choice of the nodes {ci}, allows to approximate the given integral to any degree 
of accuracy, even though, when using finite precision arithmetic, it suffices to approximate it to 
machine precision. We observe that, since the quadrature is exact for polynomials of degree q — 1, 
then its remainder depends on the q-th derivative of the integrand with respect to r. Consequently, 
considering that P^ l \c) = 0, for i > j, one has 

1 & 
A 3 -(fc) = / P J (T)f(u(to + Th))dT-J2b e P ] (c e )f(u(to + c e h)) = 0(h q - J ), (3.3) 
Jo e=1 

j = 0, . . . , r - 1. Thus, (|3~Tj) is equivalent to the ODE, 

r-1 „1 

u'(t +ch)=y2P j (c)('y j {u)-A j (h)), ce[0,l], 7j -(«)= / Pj{T)f{u{t Q +Th))dr, (3.4) 

with u(to) = yo, in place of (|2.3[) . The following result then holds true. 

Theorem 3.1. In the above hypotheses: y(to + h) — u(to + h) = 0(h p+1 ), with p — min(g,2r). 

Proof. The proof is quite similar to that of Theorem l2.2l by virtue of Lemma l2.1l and (|3.3[) - (|3.4I) . 
one obtains 

y(t + h) - u(t + h) = y(t + h,t ,y ) -y{t + h,t + h,u(t + h)) 

r to+h d rt +h / q q \ 

-y(t + h,r,u(r))dT = / — y(t + h,r, u(r)) + — y(t + /i,r,u(r))u'(r) dr 



o dr J to \dr du 

= h ( $(t + h,t + ch, u{t + ch)) (-f(u(t + ch)) + u'(t + ch)) dc 
Jo 

„1 /r-1 oo \ 

= h / <f>(t + h,t + ch,u(t +ch)) V Pj-(c) Aj-(/i) -V 7j(u)Pj(c) dc 
Jo \ 3=0 j=r J 

= ( [ Pj(T)$(t + h,t + ch, u(t + ch)) dc) Aj(u) 

3=0 ' 

(/ ^'( T ) $ (*° + Mo + e/i, it(*o + c/i)) dc^ 7j(u) 

r— 1 oo 

= hY,0(h j )0(h q - j ) - hY,0(h j )0(h j ) = 0{h q+1 ) + 0{h 2r+1 ). 

3=0 3=r 

□ 

As an immediate consequence, one has the following result. 

Corollary 2. Lei q be the order of the quadrature formula defined by the abscissae {ci}. Then, 
the order of the method (3. 2\) for approximating (1.1]) . with y± = u(to + h), is p = min(q, 2r). 

Concerning the linear stability analysis, by considering that a quadrature formula of order q > 2r 
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is exact when the integrand is a polynomial of degree at most 2r — 1, the following result immediately 
derives from Theorem 12.31 

COROLLARY 3. Let q be the order of the quadrature formula defined by the abscissae {cj}. If 
q > 2r, then the method h3.2\) is perfectly A-stable^ 

In the case r = 1, the above results apply to the methods in [10] (see also 

3.1. Runge-Kutta formulation. By setting, as usual, Ui = u(to + Cih), fi = f(ui), i = 
1, . . . , fc, p.2p can be rewritten as 

r—1 pa k 

Ui ^Vo + h^ / Pj(r) dr ^ bePj(ce)fe, i = l,...,k. (3.5) 

j=o ^° fci 

Moreover, since q > r > degw, one has y\ = u(t + h) = y + h^^iWfi- Consequently, the 
methods which Corollary [2] refers to are the subclass of Runge-Kutta methods with the following 
tableau: 



Cl 






A = (a,,) ee ^ Efco Peicj) J C< Pt(r) dr) 


Ck 






bi ... b k 



(3.6) 



In particular, in [4] it has been proved that when the nodes {cj} coincide with the k Gauss points 
in the interval [0,1], then 

A = APV T n, 

where A € M. kxk is the matrix in the Butcher tableau of the fc-stages Gauss method, V = (Pj-i (ci)) G 
M' £Xr , and il = diag(6i, . . . , bk)- In such a way, when k = r, one obtains the classical r-stages Gauss 
collocation method. Consequently, (|3.6j) can be regarded as a generalization of the classical Runge- 
Kutta collocation methods, (|3.2[) being interpreted as extended collocation conditions. 

3.2. Hamiltonian Boundary Value Methods (HBVMs). When considering a canonical 
Hamiltonian problem (|1.2|) . the discretization of the integrals appearing in (|2.3[) by means of a 
Gaussian formula at k nodes results in the HBVM(/c, r) methods introduced in [3] El For such 
methods we derive, in a novel way with respect to [TJ [2j [3], the following result. 

COROLLARY 4. For all k >r, HBVM(k,r) is perfectly A-stable and has order 2r. The method 
is energy conserving for all polynomial Hamiltonians of degree not larger than 2k/r. 



I.e., its absolute stability region coincides with the left- half complex plane, C , [§]. 
2 A different discretization, based at k + 1 Lobatto abscissae, was previously considered in [2]. 
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Proof. The result on the order and linear stability easily follow from Corollaries [5] and [3J 
respectively. Concerning energy conservation, one has 



H{u{t + h)) - H(y ) 



VH(u(t + Th)) T u'(t + rh) dr 

k 



*' j=o 1=1 



c t h)) dr 



J=0 



provided that 



=1 

i T 



Pj(T)VH(u(t +Th))dT 



J 



J2 b ePj(ce)VH{u{t + c e h)) 



U=l 



„1 k 

/ Pj{T)\7H{u{t + rh)) dr = V&^(Q)Wf Wt + cth)). 
Jo , , 



= 0, 



(3.7) 



In the case where H is a polynomial of degree v, this is true provided that the integrand is a 
polynomial of degree at most 2k — 1. Consequently, i>r — 1 < 2k — 1, i.e., v < 2k/ r. □ 

In the case of general Hamiltonian problems, if we consider the limit as k — > oo we recover 
formulae (|2.3[) . which have been called HBVM(oo,r) (or, more in general, oo-HBVMs) [IJ[3]: by the 
way, (|2.4[) is nothing but the Master Functional Equation in [TJ|3]. In the particular case r = 1, we 
derive the averaged vector field in [14j (for a related approach see [7]). 

Remark 2. We observe that, in the case of polynomial Hamiltonian systems, if \3. 7[ ) holds true 
for k = k* , then 

HBVM(fc,r) = HBVM(F,r) = HBVM(oo, r), Vfc > k* . 



That is, 113.1]) coincides with V2.S\) , for all k > k* . In the non polynomial case, the previous con- 
clusions continue "practically" to hold, provided that the integrals are approximated within machine 
precision. 

Remark 3. As is easily deduced from the arguments in Section \3.1[ the HBVM(r,r) method is 
nothing but the r -stages Gauss method of order 2r (see also 

4. Numerical Tests. We here provide a few numerical tests, showing the effectiveness of 
HBVMs, namely of the methods obtained in the new framework, when the problem (jl.lj) is in the 
form (|1.2j) . In particular, we consider the Kepler problem, whose Hamiltonian is (see, e.g., jH p. 9]): 

H ([ Ql , q 2 , Pl , P2 ] T ) = \ (pi +pI) {q{ + g^ . 

When started at 

(l-e, 0, 0, Va + eVa-e) ) T , ee[0,l), 

it has an elliptic periodic orbit of period 2-k and eccentricity e. When e is close to 0, the problem 
is efficiently solved by using a constant stepsize. However, it becomes more and more difficult as 
e — > 1, so that a variable-step integration would be more appropriate in this case. We now compare 
the following 6-th order methods for solving such problem over a 1000 periods interval: 

• HBVM(3,3), i.e., the GAUSS6 method, which is a symmetric and symplectic method; 
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Fig. 4.2. Kepler problem, e = 0.99, Hamiltonian (left plot) and solution (right plot) errors over 1000 periods with 
a variable stepsize, Tol = 10~ 10 . Note that HBVM(4,3) is not energy preserving, on the contrary of HBVM(15,3). 



• HBVM(4,3), which is symmetric [2] but not symplcctic nor energy preserving, since the 
Gauss quadrature formula of order 8 is not enough accurate, for this problem; 

• HBVM(15,3), which is practically energy preserving, since the Gauss formula of order 30 is 
accurate to machine precision, for this problem. 

In the two plots in Figure |4~T1 we see the obtained results when e = 0.6 and a constant stepsize is 
used: as one can see, the Hamiltonian is approximately conserved for the GAUSS6 and HBVM(4,3) 
methods, and (practically) exactly conserved for the HBVM(15,3) method. Moreover, all methods 
exhibit the same order (i.e., 6), with the error constant of the HBVM(4,3) and HBVM(15,3) methods 
much smaller than that of the symplectic GAUSS6 method. On the other hand, when e = 0.99, we 
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consider a variable stepsize implementation with the following standard mesh-selection strategy: 



( 



Tol 



) 



i/(p+i) 



h 



'new 



= 0.7-h. 



n 



err n 



where p = 6 is the order of the method, Tol is the used tolerance, h n is the current stepsize, and 
err n is an estimate of the local error. According to what stated in the literature, see, e.g., |15[ 
pp. 125-127], Q21 p. 235], [HI pp. 303-305], this is not an advisable choice for symplectic methods, 
for which a drift in the Hamiltonian appears, and a quadratic error growth is experienced, as is 
confirmed by the plots in Figure |4~21 The same happens for the method HBVM(4,3), which is not 
energy preserving. Conversely, for the (practically) energy preserving method HBVM(15,3), no drift 
in the Hamiltonian occurs and a linear error growth is observed. 

Remark 4. We observe that more refined (though more involved ) mesh selection strategies exist 
for symplectic methods aimed to avoid the numerical drift in the Hamiltonian and the quadratic error 
growth (see, e.g., J3[ Chapter VIII]). However, we want to emphasize that they are no more necessary, 
when using energy preserving methods, since obviously no drift can occur, in such a case. 

5. Conclusions. In this paper, we have presented a general framework for the derivation and 
analysis of effective one-step methods, which is based on a local Fourier expansion of the problem 
at hand. In particular, when the chosen basis is a polynomial one, we obtain a large subclass of 
Runge-Kutta methods, which can be regarded as a generalization of Gauss collocation methods. 

When dealing with canonical Hamiltonian problems, the methods coincide with the recently 
introduced class of energy preserving methods named HBVMs. A few numerical tests seem to show 
that such methods are less sensitive to a wider class of perturbations, with respect to symplectic, 
or symmetric but non energy conserving, methods. As matter of fact, on the contrary of the latter 
methods, they can be conveniently associated with standard mesh selection strategies. 
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